package scu.maqiang.homogenization;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L21;
import scu.maqiang.fes.FES1L31;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DoubleIntegerPair;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class SOTS1DL21 {

	public SOTS1DL21(FES1L21 fsL2, FES1L31 fsL3) {
		this.fsL2 = fsL2;
		this.fsL3 = fsL3;
	}
	
	/**
	 * һά�ȴ������ⵥ����������Ȼ�ϵ�����
	 * @param kxxFunc �ȴ���ϵ������
	 * @param param ����
	 * @param FOCS �õ�һ�׵�������
	 * @param D �õ����Ȼ��ȴ���ϵ��
	 * @param SOCS �õ����׵�������
	 */
	public void KxxCellSolution(ScalarFunc[] kxxFunc, double[] param, double[][] FOCS, double[][] D,
			double[][] SOCS) {
		SRMatrix A = new SRMatrix(fsL2.GetNdof());
		fsL2.assembleStiff(kxxFunc, param, BVPType.COMMON, A);
		double[] RHS = new double[fsL2.GetNdof()];
		fsL2.assembleGlobalVector(kxxFunc, param, SOTS1DL21::elementkxxFirstOrder, BVPType.COMMON, RHS);
		fsL2.applyBC_MBN(A, Direct.All, 1, 2);
		fsL2.applyBC_MBN(RHS, Direct.All, 0.0, 1, 2);
		NewIterSSolver solver = new NewIterSSolver(A);
		MVO.fill(FOCS, 0.0);
		solver.PCGSSOR(RHS, FOCS[0], 1.5, 1);
//
		computeHomogenizedKxx(kxxFunc, param, FOCS, D);
//
		double[][] RHS2 = new double[1][fsL2.GetNdof()];
//		fsT3.assembleGlobalVector(kxxFunc, param, FOCS, D, SOTS2DT361::elementKxxSecondOrder, BVPType.COMMON, RHS);

//		public void assembleGlobalVector(ScalarCoefFunc[] coeffFunc, double[] param, double[][] coef1, double[][] coef2,
//				ElementMMVectorFunc eleFunc, BVPType tp, double[][] RHS)
		fsL2.assembleGlobalVector(kxxFunc, param, FOCS, D, SOTS1DL21::elementKxxSecondOrder, BVPType.COMMON, RHS2);
		MVO.fill(SOCS, 0.0);
		fsL2.applyBC_MBN(RHS2[0], Direct.All, 0.0, 1, 2);
		solver.PCGSSOR(RHS2[0], SOCS[0], 1.5, 1);
	}

	public static void elementkxxFirstOrder(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[] eleXY) {
		double meanF = MVO.average(coef[0]);
		eleXY[0] = meanF;
		eleXY[1] = -meanF;
	}
	
	public void computeHomogenizedKxx(ScalarFunc[] kxxFunc, double[] param, double[][] FOCS, double[][] D) {
		double[] dx = new double[2];
		double[][] xy = new double[2][1];
		Mesh1L2 mesh = (Mesh1L2) fsL2.getMesh();
		double[] NN1 = new double[2];
		double N1Dx;
		double allLength = 0.0;
		double mKxx = 0.0;
		double[][] coef = new double[kxxFunc.length][2];
		double length = 0.0;
		MVO.fill(D, 0.0);
		for (int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			int label = mesh.getElementLabel(i);
			mesh.getCoordPerElement(ele, xy);
			MVO.getSlice(FOCS[0], ele, NN1);
			length = Math.abs(xy[0][0] - xy[1][0]);
			dx[0] = -1.0 / length;
			dx[1] = 1.0 / length;
			N1Dx = MVO.dot_product(NN1, dx);
			allLength += length;
			fsL2.getCoefFromFunc(kxxFunc, xy, label, param, coef);
			mKxx = MVO.average(coef[0]);
			D[0][0] += mKxx * (1 + N1Dx) * length;
		}
		D[0][0]  /= allLength;
	}
	
	public static void elementKxxSecondOrder(int[] ele, double[][] coord, double[][] kxx, double[][] FOCS, double[][] homo, BVPType tp,
			double[][] eleXY) {
		double[] x = MVO.col(coord, 0);
		double L = Math.abs(x[1] - x[0]);
		double[] dx = {-1.0 /L, 1.0 / L};
		double N1Dx = MVO.dot_product(dx, FOCS[0]);
		double mKxx = MVO.average(kxx[0]);
		for (int i = 0; i < 2; i++) {
			eleXY[0][i] = (mKxx - homo[0][0] + mKxx * N1Dx - mKxx * MVO.sum(FOCS[0]) * dx[i]) * L / 2.0; // 11
		}
	}
	
	public void assembleHeatSOTSSolution(FES1L21 fsFine, double[] T_0, double[] T_1, double[] T_2, FES1L21 fsHomo,
                                         double[] T_Homo, FES1L21 fsCell, double[][] FOCS, double[][] SOCS, int N, double ep) {
		double[] xy = new double[1];
		double[] cellXY = new double[1];
		
		double[] eleValHomo = new double[2];
		double[] eleValCell = new double[2];

		double[][] grad = new double[1][fsHomo.GetNdof()];
		double[][] Se = new double[1][fsHomo.GetNdof()];
		//fsHomo.computeGradient(T_Homo, grad);
		
		//fsHomo.computeFirstAndSecondDerivative(T_Homo, grad, Se);
		fsHomo.computeGradient(T_Homo, grad);
		fsHomo.computeGradient(grad[0], Se);

		double corrector1, corrector2;
		
		DoubleIntegerPair homoData = new DoubleIntegerPair(0.0, 0);
		int[] eleIdxHomo = fsHomo.mesh.getElement(0).clone();
		
		DoubleIntegerPair cellData = new DoubleIntegerPair(0.0, 0);
		int[] eleIdxCell = fsCell.mesh.getElement(0).clone();
		double val = 0.0;
		for (int i = 0, n = fsFine.GetNdof(); i < n; i++) {
			fsFine.mesh.getCoordPerNode(i, xy);
			cellXY[0] = xy[0] * N - (int) (xy[0] * N);
			fsHomo.getValue(T_Homo, xy, homoData, eleIdxHomo, eleValHomo);
			T_0[i] = homoData.first;
			fsCell.getValue(FOCS[0], cellXY, cellData, eleIdxCell, eleValCell);
			val = cellData.first;

			corrector1 = val * MVO.dot_product(grad[0], eleIdxHomo, eleValHomo);
			T_1[i] = T_0[i] + ep * corrector1;
			
			corrector2 = MVO.dot_product(SOCS[0], eleIdxCell, eleValCell) * MVO.dot_product(Se[0], eleIdxHomo, eleValHomo);
			T_2[i] = T_1[i] + ep * ep * corrector2;
		}
	}
	

	
	FES1L21 fsL2;
	FES1L31 fsL3;
}
